Baby Mental Life: Study 2

Exploratory analyses

“Baby Mental Life: Study 2” was conducted on MTurk on 2018-08-04.

Our planned sample was 300 participants, and we anticipated that roughly 80% of recruited participants would pass all of our attention checks, so we initially recruited 378 participants (on the idea that ~80% of 378 ~ 300 participants; note that for administrative purposes we need to recuit participants in batches that were divisible by 9). After filtering out participants who failed at least one of our attention checks, we ended up retaining fewer than 300 participants, so we recruited an additional 16 participants for a total of 394 people recruited. At each stage, we recruited women and men through separate studies, in hopes of acquiring a roughly equal split between genders.

In the end, we ended up with a sample of 304 participants who passed our attention checks, 237 of whom came from unique GPS coordinates.

For this first pass, these data exclude participants where there is another participant with an identical set of GPS coordinates as recorded by Qualtrics.

Each participant assessed children’s mental capacities at 13 target ages between the ages of 0 and 5 years. For each target, they rated 20 mental capacities on a scale from 0 (not at all capable) to 100 (completely capable).

For more details about the study, see our preregistration here.

Here we run some exploratory analyses on these data.

# load required libraries
library(tidyverse)
library(langcog) # source: https://github.com/langcog/langcog-package
library(psych)
library(lme4)
# set theme for ggplots
theme_set(theme_bw())

# run source code (extra home-made functions)
source("./scripts/max_factors_efa.R")
source("./scripts/plot_fun.R")
source("./scripts/reten_fun.R")
source("./scripts/table_fun.R")
source("./scripts/data_prep.R")

NAs introduced by coercionattributes are not identical across measure variables;
they will be droppedJoining, by = "question_qualtrics"

Treating factors as categories

Preliminaries

First, I’ll do the factor analysis, and get a list of the top 5 items by factor:

# load in S1 efa in case we need it
efa_S1 <- readRDS("../study 1/s1_efa.rds")
demo_S1 <- read.csv("../study 1/s1_demo.csv")

# conduct S2 efa
efa_S2 <- fa(d_all, nfactors = 4, rotate = "oblimin", fm = "minres",
             scores = "tenBerge", impute = "median")

Loading required namespace: GPArotation

table_fun(efa_S2, num_items = 5, pos_abs = "abs")

capacity	loading
Factor 1
feeling_distressed	0.82
feeling_overwhelmed	0.82
feeling_frustrated	0.78
feeling_helpless	0.76
feeling_lonely	0.60
Factor 2
having_self_control	0.96
controlling_their_emotions	0.93
telling_right_from_wrong	0.93
planning	0.89
reasoning_about_things	0.89
Factor 3
getting_hungry	0.83
feeling_pain	0.79
feeling_tired	0.67
hearing_sounds	0.58
feeling_physically_uncomfortable	0.49
Factor 4
feeling_happy	0.82
finding_something_funny	0.81
feeling_excited	0.79
loving_somebody	0.64
learning_from_other_people	0.51

Now I’ll use this to define categories of mental capacities:

# get items by factor
factors_S2 <- efa_S2$loadings[] %>%
  data.frame() %>%
  rownames_to_column("capacity") %>%
  gather(factor, loading, -capacity) %>%
  group_by(capacity) %>%
  top_n(1, loading) %>%
  ungroup() %>%
  select(-loading) %>%
  mutate(factor_names = recode_factor(factor,
                                      "MR1" = "Negative emotions",
                                      "MR2" = "Cognition & control",
                                      "MR3" = "Bodily sensations",
                                      "MR4" = "Positive/social emotions"),
         factor = factor(factor))
# make new dataframe
d_cat <- d_all %>%
  rownames_to_column("subid_target") %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(capacity, response, -c(subid, starts_with("target"))) %>%
  left_join(factors_S2) %>%
  left_join(d_demo %>%
              select(ResponseId, Parent) %>%
              rename(subid = ResponseId, parent = Parent) %>%
              mutate(subid = as.character(subid)))

Joining, by = "capacity"
Joining, by = "subid"

And get an average “score” for each of these factors for each participant:

d_cat_scored <- d_cat %>%
  group_by(subid, parent, 
           target, target_num, target_ord, 
           factor, factor_names) %>%
  summarise(score = mean(response, na.rm = T)) %>%
  ungroup()

Plots

# bootstrapped means and 95% CIs
d_cat_scored_boot <- d_cat_scored %>%
  group_by(target, target_num, target_ord, factor, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup()

ggplot(d_cat_scored,
       aes(x = target_num, y = score, color = factor_names)) +
  facet_grid(~ factor_names) +
  geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_scored_boot, color = "black",
            aes(y = mean, group = factor_names)) +
  geom_pointrange(data = d_cat_scored_boot, color = "black", fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Set2", guide = "none") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, untransformed",
       x = "target age (months)")

ggplot(d_cat_scored,
       aes(x = target_num, y = score, color = factor_names)) +
  facet_grid(~ factor_names) +
  geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_scored_boot, color = "black",
            aes(y = mean, group = factor_names)) +
  geom_pointrange(data = d_cat_scored_boot, color = "black", fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Set2", guide = "none") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, square-root-transformed",
       x = "age after square-root transformation (months)")

ggplot(d_cat_scored,
       aes(x = target_ord, y = score, color = factor_names)) +
  facet_grid(~ factor_names) +
  geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_scored_boot, color = "black",
            aes(y = mean, group = factor_names)) +
  geom_pointrange(data = d_cat_scored_boot, color = "black", fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Set2", guide = "none") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Age as an ordinal variable",
       x = "age (ordinal)")

Regressions

# contrasts(d_cat_scored$target_ord) <- contr.poly(13)
# contrasts(d_cat_scored$factor) <- contr.sum(4)
contrasts(d_cat$target_ord) <- contr.poly(13)
contrasts(d_cat$factor) <- contr.sum(4)

r1 <- lmer(response ~ factor * poly(target_num, 3) +
             (1 + factor + target_num | subid),
           d_cat %>%
             mutate(target_num = scale(target_num, center = F)))
summary(r1)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ factor * poly(target_num, 3) + (1 + factor + target_num |  
    subid)
   Data: d_cat %>% mutate(target_num = scale(target_num, center = F))

REML criterion at convergence: 536724.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.5360 -0.4298  0.0225  0.5167  5.5499 

Random effects:
 Groups   Name        Variance Std.Dev. Corr                   
 subid    (Intercept) 200.49   14.159                          
          factor1     156.13   12.495    0.65                  
          factor2     145.83   12.076   -0.44 -0.59            
          factor3      69.75    8.352   -0.55 -0.34 -0.24      
          target_num   51.83    7.199   -0.64 -0.42  0.50  0.31
 Residual             334.01   18.276                          
Number of obs: 61620, groups:  subid, 237

Fixed effects:
                               Estimate Std. Error t value
(Intercept)                     68.3739     0.7634  89.569
factor1                          6.9512     0.8216   8.461
factor2                        -38.7255     0.7947 -48.728
factor3                         25.1428     0.5573  45.116
poly(target_num, 3)1          3241.2664    88.4095  36.662
poly(target_num, 3)2         -1237.5812    18.2759 -67.717
poly(target_num, 3)3           532.0777    18.2759  29.114
factor1:poly(target_num, 3)1  -544.0136    31.6548 -17.186
factor2:poly(target_num, 3)1  2670.8981    31.6548  84.376
factor3:poly(target_num, 3)1 -2550.8843    31.6548 -80.584
factor1:poly(target_num, 3)2    32.5828    31.6548   1.029
factor2:poly(target_num, 3)2    38.2604    31.6548   1.209
factor3:poly(target_num, 3)2   881.6516    31.6548  27.852
factor1:poly(target_num, 3)3    48.0128    31.6548   1.517
factor2:poly(target_num, 3)3  -391.9779    31.6548 -12.383
factor3:poly(target_num, 3)3  -354.1489    31.6548 -11.188

Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE)  or
     vcov(x)     if you need it

This is weirdly hard to model… lots of alternative models didn’t converge. Also, those extremely high t-values make me a little suspicious, as do the identical standard errors for the last 9 terms. So, let’s take this with a bit of a grain of salt for now.

Parents vs. non-parents

Study 2

Preliminaries

d_parent <- d_all %>%
  rownames_to_column("subid_target") %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(capacity, response, -c(subid, starts_with("target"))) %>%
  left_join(d_demo %>%
              select(ResponseId, Parent) %>%
              rename(subid = ResponseId, parent = Parent) %>%
              mutate(subid = as.character(subid)))

Joining, by = "subid"

parent_counts <- d_demo %>% 
  distinct(ResponseId, Parent) %>% 
  rename(subid = ResponseId, parent = Parent) %>%
  count(parent) %>% 
  mutate(proportion = n/sum(n))
knitr::kable(parent_counts, digits = 3)

parent	n	proportion
No	146	0.616
Yes	88	0.371
NA	3	0.013

d_demo %>%
  count(Parent, ChildrenOldestAge_collapse)

Given these numbers, I think our best bet it just to look at parents vs. non-parents, and not try to separate out parents of young children (too few!).

scores_S2 <- efa_S2$scores %>%
  data.frame() %>%
  rownames_to_column("subid_target") %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(factor, score, -c(subid, starts_with("target"))) %>%
  left_join(d_parent %>% distinct(subid, parent)) %>%
  left_join(d_cat %>% distinct(factor, factor_names)) %>%
  mutate(factor = factor(factor))

Joining, by = "subid"
Joining, by = "factor"
Column `factor` joining character vector and factor, coercing into character vector

Plots

Factor scores

# bootstrapped means and 95% CIs
d_parent_boot <- scores_S2 %>%
  filter(!is.na(parent)) %>%
  group_by(parent, target, target_num, target_ord, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup() %>%
  data.frame()

scores_S2 %>%
  filter(!is.na(parent)) %>%
  ggplot(aes(x = target_num, y = score, color = parent)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.05) +
  geom_line(data = d_parent_boot, aes(y = mean, group = parent)) +
  geom_pointrange(data = d_parent_boot, fatten = 1.5, # note: too close to dodge
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = parent), #color = "black",
  #             method = "lm", formula = "y ~ poly(x, 3)") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
  labs(title = "Developmental trajectories of factor scores, by parent status",
       subtitle = "Exact age, untransformed",
       color = "parent status: ",
       x = "target age (months)", y = "score (NOTE: scales differ)")

scores_S2 %>%
  filter(!is.na(parent)) %>%
  ggplot(aes(x = target_num, y = score, color = parent)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.05) +
  geom_line(data = d_parent_boot, aes(y = mean, group = parent),
            position = position_dodge(width = 0.3)) +
  geom_pointrange(data = d_parent_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper),
                  position = position_dodge(width = 0.3)) +
  # geom_smooth(aes(group = parent), #color = "black",
  #             method = "lm", formula = "y ~ poly(x, 3)") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
  scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
  labs(title = "Developmental trajectories of factor scores, by parent status",
       subtitle = "Exact age, square-root transformed",
       color = "parent status: ",
       x = "age after square-root transformation (months)", 
       y = "score (NOTE: scales differ)")

scores_S2 %>%
  filter(!is.na(parent)) %>%
  ggplot(aes(x = target_ord, y = score, color = parent)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.05) +
  geom_line(data = d_parent_boot, aes(y = mean, group = parent),
            position = position_dodge(width = 0.5)) +
  geom_pointrange(data = d_parent_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper),
                  position = position_dodge(width = 0.5)) +
  # geom_smooth(aes(group = parent), #color = "black",
  #             method = "lm", formula = "y ~ poly(x, 3)") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
  labs(title = "Developmental trajectories of factor scores, by parent status",
       subtitle = "Age as ordinal variable",
       color = "parent status: ",
       x = "age (ordinal)", 
       y = "score (NOTE: scales differ)")

Category ‘summary scores’

# bootstrapped means and 95% CIs
d_cat_parent_scored_boot <- d_cat_scored %>%
  filter(!is.na(parent)) %>%
  group_by(parent, target, target_num, target_ord, factor, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup()

ggplot(d_cat_parent_scored_boot,
       aes(x = target_num, y = score, color = parent)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_parent_scored_boot, 
            aes(y = mean, group = parent)) +
  geom_pointrange(data = d_cat_parent_scored_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, untransformed",
       x = "target age (months)",
       y = "score (NOTE: scales differ)")

ggplot(d_cat_parent_scored_boot,
       aes(x = target_num, y = score, color = parent)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_parent_scored_boot, 
            aes(y = mean, group = parent)) +
  geom_pointrange(data = d_cat_parent_scored_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, square-root transformation",
       x = "age after square-root transformation (months)",
       y = "score (NOTE: scales differ)")

ggplot(d_cat_parent_scored_boot,
       aes(x = target_ord, y = score, color = parent)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_parent_scored_boot, 
            aes(y = mean, group = parent)) +
  geom_pointrange(data = d_cat_parent_scored_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Age as an ordinal variable",
       x = "age (ordinal)",
       y = "score (NOTE: scales differ)")

All of these visualizations suggest that, relative to non-parents (n = NA), parents (n = NA) tended to perceive greater abilities, in all domains except for bodily sensations. For both negative emotions and positive/social emotions, this difference is greatest in the mid-range of the target ages (12-24 months); for cognition & control, it’s greater at the older end of the target ages (e.g., 3-5 years). It doesn’t look like a huge effect, but it’s intruiging.

Regressions

There are many ways that we could choose to model this - I’ll try out the factor scores first, adapting our primary regression models (not in this notebook).

contrasts(scores_S2$parent) <- contr.treatment(2, base = 1) # baseline: NON-parents
contrasts(scores_S2$factor) <- contr.sum(4)

r2 <- lmer(score ~ target_num * factor * parent
           + (target_num + factor | subid),
           scores_S2 %>%
             mutate(target_num = target_num/12))
summary(r2)

Linear mixed model fit by REML ['lmerMod']
Formula: score ~ target_num * factor * parent + (target_num + factor |  
    subid)
   Data: scores_S2 %>% mutate(target_num = target_num/12)

REML criterion at convergence: 22480.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-9.3553 -0.4388  0.0524  0.5114  5.2974 

Random effects:
 Groups   Name        Variance Std.Dev. Corr                   
 subid    (Intercept) 0.3491   0.5908                          
          target_num  0.0165   0.1284   -0.73                  
          factor1     0.1947   0.4413    0.49 -0.40            
          factor2     0.2303   0.4799   -0.66  0.56 -0.40      
          factor3     0.2605   0.5104    0.24  0.00 -0.22 -0.61
 Residual             0.2919   0.5403                          
Number of obs: 12168, groups:  subid, 234

Fixed effects:
                            Estimate Std. Error t value
(Intercept)                -0.424399   0.049602  -8.556
target_num                  0.279097   0.011334  24.624
factor1                     0.078297   0.039260   1.994
factor2                    -0.273518   0.042253  -6.473
factor3                     0.253010   0.044634   5.669
parent2                     0.070601   0.080884   0.873
target_num:factor1         -0.058596   0.006811  -8.603
target_num:factor2          0.190717   0.006811  28.001
target_num:factor3         -0.170985   0.006811 -25.104
target_num:parent2          0.003333   0.018482   0.180
factor1:parent2             0.005619   0.064020   0.088
factor2:parent2            -0.051089   0.068900  -0.741
factor3:parent2             0.031692   0.072783   0.435
target_num:factor1:parent2 -0.002939   0.011107  -0.265
target_num:factor2:parent2  0.053753   0.011107   4.840
target_num:factor3:parent2 -0.046204   0.011107  -4.160

Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE)  or
     vcov(x)     if you need it

The last three coefficients here are where the action is: According to this model, the difference between parents and non-parents in perceived developmental trajectories (target_num:parent2, which is not significant collapsing across factors) is exaggerated in the domain of cognition and control (target_num:factor2:parent2) and diminished in the domain of bodily sensations (target_num:factor3:parent2).

Let’s try adding polynomial effects:

r3 <- lmer(score ~ poly(target_num, 3) * factor * parent
           + (poly(target_num, 1) + factor | subid),
           scores_S2 %>%
             mutate(target_num = target_num/12))
summary(r3)

Linear mixed model fit by REML ['lmerMod']
Formula: score ~ poly(target_num, 3) * factor * parent + (poly(target_num,  
    1) + factor | subid)
   Data: scores_S2 %>% mutate(target_num = target_num/12)

REML criterion at convergence: 19808.5

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-10.3042  -0.4652   0.0017   0.5305   5.8317 

Random effects:
 Groups   Name                Variance Std.Dev. Corr                   
 subid    (Intercept)           0.2278  0.4773                         
          poly(target_num, 1) 519.8012 22.7991  -0.51                  
          factor1               0.1982  0.4452   0.45 -0.39            
          factor2               0.2338  0.4835  -0.60  0.54 -0.40      
          factor3               0.2640  0.5138   0.30  0.00 -0.22 -0.61
 Residual                       0.2316  0.4812                         
Number of obs: 12168, groups:  subid, 234

Fixed effects:
                                       Estimate Std. Error t value
(Intercept)                           -0.030561   0.039884  -0.766
poly(target_num, 3)1                  48.855288   1.983984  24.625
poly(target_num, 3)2                 -19.550418   0.613116 -31.887
poly(target_num, 3)3                   8.448995   0.613116  13.780
factor1                               -0.004388   0.038065  -0.115
factor2                               -0.004395   0.041145  -0.107
factor3                                0.011731   0.043587   0.269
parent2                                0.075304   0.065038   1.158
poly(target_num, 3)1:factor1         -10.257028   1.061949  -9.659
poly(target_num, 3)2:factor1           3.413024   1.061949   3.214
poly(target_num, 3)3:factor1          -1.785482   1.061949  -1.681
poly(target_num, 3)1:factor2          33.384462   1.061949  31.437
poly(target_num, 3)2:factor2           3.384910   1.061949   3.187
poly(target_num, 3)3:factor2          -6.433547   1.061949  -6.058
poly(target_num, 3)1:factor3         -29.930494   1.061949 -28.184
poly(target_num, 3)2:factor3           8.594647   1.061949   8.093
poly(target_num, 3)3:factor3          -3.251241   1.061949  -3.062
poly(target_num, 3)1:parent2           0.583455   3.235228   0.180
poly(target_num, 3)2:parent2          -2.031708   0.999792  -2.032
poly(target_num, 3)3:parent2           1.883539   0.999792   1.884
factor1:parent2                        0.001472   0.062072   0.024
factor2:parent2                        0.024763   0.067094   0.369
factor3:parent2                       -0.033507   0.071076  -0.471
poly(target_num, 3)1:factor1:parent2  -0.514473   1.731690  -0.297
poly(target_num, 3)2:factor1:parent2  -0.041808   1.731690  -0.024
poly(target_num, 3)3:factor1:parent2   1.099314   1.731690   0.635
poly(target_num, 3)1:factor2:parent2   9.409337   1.731690   5.434
poly(target_num, 3)2:factor2:parent2  -2.361622   1.731690  -1.364
poly(target_num, 3)3:factor2:parent2  -1.057673   1.731690  -0.611
poly(target_num, 3)1:factor3:parent2  -8.087893   1.731690  -4.671
poly(target_num, 3)2:factor3:parent2   6.669873   1.731690   3.852
poly(target_num, 3)3:factor3:parent2  -4.830063   1.731690  -2.789

Correlation matrix not shown by default, as p = 32 > 12.
Use print(x, correlation=TRUE)  or
     vcov(x)     if you need it

A lot to sort through here, but here are some observations:

• Collapsing across factors, it seems like parents might perceive exaggerated non-linearities? (poly(target_num, 3)2:parent2 and poly(target_num, 3)3:parent2)
• Compared to collapsing across factors, not a lot of additional parent vs. non-parent differences in negative emotions (factor 1)
• Compared to collapsing across factors, parents perceived more dramatic (linear) growth in cognition/control (factor 2; poly(target_num, 3)1:factor2:parent2)
• Compared to collapsing across factors, in the domain of bodily sensations (factor 3), parents perceived less growth (poly(target_num, 3)1:factor3:parent2). I’m not sure how to interpret the differences in non-linearity here.

You could imagine re-running any of these with a square-root transformation on target age, and/or with the ‘summary scores’ by category, but I’m not going to do that right now.

Study 1

Preliminaries

parent_counts_S1 <- demo_S1 %>% 
  distinct(ResponseId, Parent) %>% 
  rename(subid = ResponseId, parent = Parent) %>%
  count(parent) %>% 
  mutate(proportion = n/sum(n))
knitr::kable(parent_counts_S1, digits = 3)

parent	n	proportion
No	147	0.653
Yes	77	0.342
NA	1	0.004

scores_S1 <- efa_S1$scores %>%
  data.frame() %>%
  rownames_to_column("subid_target") %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(factor, score, -c(subid, starts_with("target"))) %>%
  left_join(demo_S1 %>% distinct(ResponseId, Parent) %>%
              rename(subid = ResponseId, parent = Parent) %>%
              mutate(subid = as.character(subid))) %>%
  mutate(factor_names = recode_factor(factor,
                                      "MR4" = "Negative emotions (S1 F4)",
                                      "MR1" = "Cognition & control (S1 F1)",
                                      "MR2" = "Bodily sensations (S1 F2)",
                                      "MR3" = "Positive/social emotions (S1 F3)"),
         factor = factor(factor))

Joining, by = "subid"

Plots

Factor scores

# bootstrapped means and 95% CIs
d_parent_boot_S1 <- scores_S1 %>%
  filter(!is.na(parent)) %>%
  group_by(parent, target, target_num, target_ord, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup() %>%
  data.frame()

scores_S1 %>%
  filter(!is.na(parent)) %>%
  ggplot(aes(x = target_num, y = score, color = parent)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.05) +
  geom_line(data = d_parent_boot_S1,
            position = position_dodge(width = 5),
            aes(y = mean, group = parent)) +
  geom_pointrange(data = d_parent_boot_S1, 
                  position = position_dodge(width = 5),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = parent), #color = "black",
  #             method = "lm", formula = "y ~ poly(x, 3)") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
  labs(title = "Developmental trajectories of factor scores, by parent status (STUDY 1)",
       subtitle = "Exact age, untransformed",
       color = "parent status: ",
       x = "target age (months)", y = "score (NOTE: scales differ)")

I won’t bother to re-plot with different treatments of age, since there are only three target ages here.

Category ‘summary scores’

# make new dataframe
d_cat_S1 <- d_all_S1 %>%
  rename(subid_target = X) %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(capacity, response, -c(subid, starts_with("target"))) %>%
  left_join(factors_S2) %>%
  left_join(demo_S1 %>%
              select(ResponseId, Parent) %>%
              rename(subid = ResponseId, parent = Parent) %>%
              mutate(subid = as.character(subid))) %>%
  filter(!is.na(factor))

Joining, by = "capacity"
Joining, by = "subid"

And get an average “score” for each of these factors for each participant:

d_cat_scored_S1 <- d_cat_S1 %>%
  group_by(subid, parent, 
           target, target_num, target_ord, 
           factor, factor_names) %>%
  summarise(score = mean(response, na.rm = T)) %>%
  ungroup()

# bootstrapped means and 95% CIs
d_cat_parent_scored_boot_S1 <- d_cat_scored_S1 %>%
  filter(!is.na(parent)) %>%
  group_by(parent, target, target_num, target_ord, factor, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup()

ggplot(d_cat_parent_scored_boot_S1,
       aes(x = target_num, y = score, color = parent)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_parent_scored_boot_S1,
            position = position_dodge(width = 5),
            aes(y = mean, group = parent)) +
  geom_pointrange(data = d_cat_parent_scored_boot_S1,
                  position = position_dodge(width = 5),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, untransformed",
       x = "target age (months)",
       y = "score (NOTE: scales differ)")

Looking quickly at these two visualizations of Study 1, I’d say that there’s a hint of the same patterns in the domains of negative emotions and positive social emotions (parents attributed more than non-parents, especially at 9 months), and in cognition/control (parents attributed more than non-parents at 5 years). So, generally consistent?

I won’t run regression analyses just now.